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Abstract 

Unsteady aerodynamic models based on windtunnel 
forced oscillation test data and analyzed with a fuzzy logic 
algorithm are incorporated into an F-16XL flight simulation 
code. The reduced frequency needed in the unsteady models 
is numerically calculated by using a limited prior time history 
of state variables in a least-square sense. Numerical examples 
are presented to show the accuracy of the calculated reduced 
frequency. Oscillatory control inputs are employed to 
demonstrate the differences in the flight characteristics based 
on unsteady and quasi-steady aerodynamic models. 
Application of the unsteady aerodynamic models is also 
presented and the results are compared with one set of F- 
16XL longitudinal maneuver flight data. It is shown that the 
main differences in dynamic response are in the lateral- 
directional characteristics, with the quasi-steady model being 
more stable than the flight vehicle, while the unsteady model 
being more unstable. Similar conclusions can also be made in 
a simulated rapid sideslipping roll. To improve unsteady 
aerodynamic modeling, it is recommended to acquire test data 
with coupled motions in pitch, roll and yaw. 

Nomenclature 

g gravitational acceleration 

H altitude 

I x ,I y ,l z - moments of inertia about the x-, y- and z-axes, 

respectively 


I XV' ^X7,' lyz 

product of inertia 

F] 1 aero, F m *ero. 

F n «xo,T,.T m ,T„ 

aerodynamic and thrust moments about 
the x-, y- and z-axes, respectively 

w 

aircraft weight 

P. q, r 

body-axis roll rate, pitch rate and yaw 
rate, respectively 

u, v, w 

aircraft velocity components along the x-, 
y-, and z-axes, respectively 

V 

total velocity 

F* IOO. Fy 

F T T T 

1 z loo* A z 

aerodynamic forces and thrust along the 
x-, y- and z-axes, respectively 

a 

angle of attack 

P 

angle of sideslip 

6., 6, 

control deflections in pitch, roll and yaw, 
respectively 

6 t 

throttle position 

<f>, 6, i)r 

Euler angles in roll, pitch, and yaw, 
respectively 

k 

reduced frequency 

b 

span length 

c 

mean aerodynamic chord 

c u 

directional cosines 


Introduction 


Aircraft dynamic characteristics have mostly been 
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studied with a quasi-steady aerodynamic model in the past 
In a quasi -stead} model, the dynamic effects on 
aerodynamic coefficients are assumed to be linearly 
proportional to the angular rates and lime rales of change 
of variables with the coefficients being typically measured 
at a fixed frequency. That is , the lag effect exhibited by the 
varying reduced frequency along a maneuvering flight path, 
and the corresponding dynamic stall effect are not properly 
accounted for. The present investigation on the unsteady 
aerodynamic effect is undertaken in view of the fact that the 
lag and dynamic stall effects may play an important role in 
aircraft stability characteristics and performance 1,2 in some 
flight conditions. 

To determine the unsteady aerodynamic effects, 
harmonic forced oscillation tests were conducted. 3 * 4 These 
test data were then analyzed to establish the aerodynamic 
models for use in the flight simulation. In the process of 
application, the most uncertain quantity to be determined is 
the reduced frequency. It is straightforward to calculate the 
reduced frequency of harmonic oscillation motion used in 
the test. However, the actual aircraft motion cannot be 
described by one single harmonic oscillation, and the 
reduced frequency of the actual motion will not be constant. 
One way to determine the reduced frequency of aircraft 
motion is an analytical method which assumes the mean 
angle of attack and the amplitude being the same as those 
used in the harmonic oscillation tests. 1 If the forced 
oscillation tests were performed at a constant mean angle of 
attack and amplitude, the analytical method should produce 
reasonable values of the “equivalent” reduced frequency. 

On the other hand, when the wind tunnel data are obtained 
at different mean angles of attack and amplitudes, this 
analytical method would not be applicable. One way to 
determine the equivalent reduced frequency is a numerical 
one. 

In the present paper, we will demonstrate one 
numerical method to determine the reduced frequency for 
use in flight simulation. The method is then employed, 
together with the fuzzy-logic aerodynamic models 
established previously 3,4 , in an F-16XL simulation code to 
investigate the unsteady aerodynamic effects on flight 
dynamics. 

Aircraft Dynamic Model 

The flight simulation is implemented by the 
following 6 -degree -of- freedom model; 


u = n>-qw+a xb 
v=pw-ru + a yb 
*=qu-pv + a lJk 
p--H x A n +H y A n +H t 
4-H z A^H y A^H, 
f=H z A^H y A n +H 2 


where 

a rJb = g( a x + C 13> 

V = s(V c b) 

a ,J, = 8(- a n + C J3> 
a x = ^F zaero )!W 

V (Jy + F y atro )IW 

a s-cr,+F^/w 


H X~P^X~^ I y~ r I a 

H r -P I x,-<l I y, +rI , 
K^Hy-qH^F^T, 
H y =pH,~rH z +F matro+ T m 
H,=rH z -qH y +F^T n 


^ = tan“ 1 (C 23 /C 33 ) 

6 - -sin -1 (C 13 ) 

V' = tan (Q2 / Cj j ) 

A* 

An=A i2 HI x J a +IJ yt )/D 

A , ,=<v,- 7 V /D 

A, =>»,>=( y* + V»> /D 

n - wv w v w 

The flight state variables (u,\\w, p,q,r) arc 
calculated by' integrating these dynamic equations. The 
following auxiliary' equations arc then used to obtain the 


(1) 


23 


33 
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angle of attack and sideslip angle and their time rates : 

a = tan' 1 (n , /w) 

& = (iAi , -H-w)/(w 2 + H' 2 ) 
p = sin M (v/F) (2 ) 

$ =(Vv- vF)/(^ 2 cosp) 

V~\ju 2 + v 2 + w 2 

Aerodynamic Models 
Quasi-steadv model 

For quasi-steady aerodynamics, the terms F XKT0 , 
«oo ) ^z»aoi f*i»ero > fmiao » ^n»ero ^ the ai re raft dynamic 
model are functions of flight state variables and control 
surface deflection angles, i.e., 

F x, y. i, «, n, at ro^ h («> P> P> ?» r, K H, 6) (3) 


where 6 stands for the control deflection angle. When the 
flight state is known, it is straightforward to obtain the 
corresponding aerodynamic forces and moments through 
interpolation. The aerodynamic effects caused by the time 
rate of angle of attack are not included in the present quasi- 
steady aerodynamic model for the F-16XL. 

Unsteady aerodynamic model 

The longitudinal and lateral -directional unsteady 
aerodynamic models were established based on the results 
of large-amplitude pitch, roll and yaw forced oscillation 
tests, and a limited amount of rotary balance data. 3,4 The 
tests were conducted in separate pitch, roll and yaw 
motions. These models were set up to have the following 
functional forms: 

Longitudinal unsteady aerodynamics model: 

CuD t m y n,Ly = / (a ,a ,a , f},S € ) (4) 

Lateral-directional unsteady aerodynamics model: 

g{a,<fi y <f>,k 2 ,y/ y y/) (5) 

where the subscripts L, I), m, n, I, y stand for lift. drag, 
pitching, yawing and rolling moments, and side force, 
respectively. The reduced frequency k, is computed from 


a(l) and q(t), representing the longitudinal unsteady 
motion characteristics On the other hand, the reduced 
frequency- k 2 is calculated from the lateral-directional 
unsteady- motion, mainly from the rolling motion if it is 
present. Note that there are no unsteady aerodynamic data 
available with the lateral-directional control surfaces 
deflected. Only the elevator effect is included in the 
longitudinal model. 

In flight simulation, the motion variables a, 0, ifr 
and their time rates are calculated at every- time instant. 

The calculation of reduced frequency kj is discussed in the 
next section. For the lateral-directional case, two reduced 
frequencies can be obtained separately from (p f , and ifr, 

y . In most cases, k 2 will be calculated from <p and ^ 
unless there is no aileron input. In this latter case, ^and 
\j/ will be used to obtain k 2 . This approach is based on the 
fact that the lateral-directional motion variables tend to 
have similar time histories. 

In addition to the aforementioned unsteady- 
aerodynamic models, the existing static aerodynamic model 
in the simulation code is used to provide the basic 
aerodynamic properties. The unsteady aerodynamic effects 
are added to the basic static aerodynamics according to the 
following expressions: 

Longitudinal unsteady aerodynamic increments: 

AC umteady\ = AC L,D.m,n,t,y = /(<*,«,«,£,, 

-f(a,a=0,a=0,k i =0,/3,S e ) {b) 

Lateral-directional unsteady aerodynamic 
increments: 

ACunsteadyl = A CL y D,m y n y t y y = g(a, <j), (f), ky , \j/ ) 

(7) 

-g(a><t>,j> = 0>k 2 =0 ,y/ y y/ = 0 ) 

Generally, <*and ^are replaced by p and r, respectively. 

Therefore, the total aerodynamic coefficients at every time 
instant are expressed by: 

^ ' L y D y m y n t £,y “ ^static unstead\ \ unsteady! ( 8 ) 

For the purpose of stability analysis, the 
derivatives of aerodynamic coefficients with respect to the 
angle of attack etc. can also be derived from the above 
equations. In addition to the numerical central difference 
method and spline interpolation to obtain these derivatives, 
the small-amplitude harmonic oscillatory method* 1 will be 
used and compared. 
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Equivalent Reduced Frequency 

In the present unsteady aerodynamic models, 
different types of forced oscillation test results (different 
mean angles of attack and different amplitudes) are used as 
source data to train the fuzzy logic models. This was done 
in the testing to accommodate a wider range of reduced 
frequency without exceeding the loading capability of the 
measuring balance. As stated previously, the angle of attack 

a, and the time rate of angle of attack & are part of the 

input parameters to calculate the unsteady aerodynamic 
coefficients in the longitudinal fuzzy logic model. The time 
history of these variables is fitted with one of a harmonic 
motion at that instant, i.e., 

a{ () -a\a cos {cot + <j > ) 

( 9 ) 

a(t) = -acosin(a)( + ^) 

where those terms on the left hand side are given and the 
unknowns are a , the local mean angle of attack, V, the 

local amplitude of the harmonic motion, <p » the phase lag, 
and co, the angular frequency. The reduced frequency k, 
and k 2 are defined as 

k x = cdclV 9 k 2 = cob 12V ( 10 ) 

The frequency <0 and other parameters are calculated 
through an optimization method by minimizing the 
following cost function (least squares) 

J - £ [a. - (a + acos(u)r. + 4>))] 2 + 

(11) 

[A. - (-aa)sin(ci)f. + <{)))] 2 

i-i 


where n is the number of dal a points used, and is assumed 
to be 20 (equivalent to 0.7 sec) in the present application. 
Those twenty points preceding the current time are 
employed in the above equation. The least-square method is 
found to converge well and gives reasonably accurate 
results. The lateral-directional equivalent reduced 
frequency is computed in the same manner. 

Numerical Results and Discussion 

Equivalent Reduced Frequency and Pilch Damping 


Assume the number of points in the motion history 
preceding the current instant is n If n > 20, it is set to 20. 
However, if n < 20, then only those available points are 
used to determine the unknowns. As is well known, the 
initial values assumed will significantly affect the 
convergence of numerical optimization. Therefore, the 
initial values of the unknowns for the current time segment 
are always taken to be the results from the previous 
segment. At the beginning, the initial values of mean a 
( a ) and amplitude ( a ) are assumed to be 35.0 deg. (the 
value used in the windtunnel testing), and the difference 
between the actual a and the assumed mean a, 
respectively. The initial values of angular frequency and 
phase angle are 1.0 and 7C, respectively. 

Figure 1 presents the comparison between the 
assumed harmonic motion with (0 = 71/2 and the predicted 
results. Except a small deviation when n < 20, the present 
method is shown to predict accurate results. The method 
has also been tested for various other assumed harmonic 
motions with similar good results (not shown). 

The convergence of the present algorithm is fast. 
For an <x(t) response with 150 points, it took only one 
second of the CPU time to calculate the reduced frequency 
at all 150 points on an Alpha VAX machine. 

The variation of longitudinal dynamic derivatives 
with the angle of attack based on the concept of small 
amplitude harmonic oscillation is presented in Figure 2. 
These dynamic derivatives are calculated by integrating the 
out-of-phase aerodynamic coefficients obtained from the 
fuzzy logic model. The amplitude is assumed to be 5 
degrees and two different reduced frequencies are presented 
in the calculation. The pitching moment is referenced to a 
pitch center at 0.558 c . The results show that predicted 
pitch damping derivatives are mostly positive because of 
the aft pitch center location. Their sign agrees with the 
hysteresis direction in that a clockwise loop produces a 
positive value for the damping derivatives. Their 
magnitude is decreased with increasing reduced frequency, 
having the same trend of variation with that calculated by 
the indicial integration method 4 5 This exercise 
demonstrates the generality and usefulness of fuzzy logic 
models as steady’ and unsteady aerody namic data base 

Dynamic Response to Oscillatory Input 

The aerodynamic moment reference point of the 
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F-16XL simulation code is al ^^=0.45^, being diflcrcnt 
from the ones used in the longitudinal and lateral- 
directional wind tunnel tests. Therefore, the pitching 
moment and yawing moment coefficients obtained from the 
fuzzy logic models must first be transferred before being 
incorporated into the simulation code by using the 
following expressions: 

~ modtl + ( X xtm2 X long^ / 1 T \ 

C — C + C (x ~ X \ f 1 

* modil ^y^ X sim2 x Ux> 

Note that x^ = 0.558 c and x tot = 0.46 c , 

To see the differences in the dynamic response 
with unsteady and quasi-steady aerodynamic models, the 
following oscillatory elevator input by the pilot is assumed: 

5 epi]ot = 15.0+25.0 sin(27T/8 t), in deg. 

The resulting dynamic response is presented in Figure 3. 
The lateral-directional state variables, such as the sideslip 
angle, roll rate and yaw rate etc., are zero. Two methods in 
calculating the dynamic derivatives with the unsteady 
aerodynamic model are presented: one with the small 
amplitude oscillation method mentioned earlier and the 
other based on a central difference method. For the latter, 
the following equation is used: 

C^ = [/(a,A,a,^,A&)-/(« -A,d,^„/?,&)]/(2A)(13) 
A = cA/(2V) 

where A= q + a if the sum is greater than 5 deg./sec.; 
otherwise, A = 5 deg./sec. For the lateral-directional 
derivatives, the following expression is used: 

Ctp = [g(a,<f>, p + Ap,k 2 ,iy,r)- p - Ap,k 2 ,y/,r)]/ (2Ap) (14) 

where Ap=3 deg./sec. and A^=bAp/(2V). The yaw 
damping derivative is calculated in a similar manner. Note 
that pitching derivatives are calculated differently to avoid 
the effect of static data when the dynamic data are not 
symmetrical about the static data. The calculated results 
indicate that the derivatives by the small amplitude method 
vary r in a more orderly manner as compared with the 
central difference method: but the order of magnitude is 
comparable. One reason for this is that in the small 
amplitude method the results depend only on a and k; 
while the central difference method depends on a as well 


in addition to a and k. One important difference between 
the unsteady model and the quasi -steady model is that the 
latter features less pitch damping at low a and higher 
damping (i.e., values being more negative) al high a, being 
opposite to what the unsteady' model indicates. The lift due 
to pitching is in general much smaller with the quasi-steady 
model. Note that during the simulation, the existing control 
system is still operational. The resulting angle of attack 
time history is seen to be the same for both the quasi-steady 
and unsteady' models. This is achieved by the control 
system with different final elevator deflection angles to 
compensate for the differences in dynamic characteristics. 
Later, a test case will be shown that this is not the case in 
other unsteady flight conditions where the control system 
designed with the quasi-steady data becomes inadequate. In 
the following, all stability derivatives will be calculated 
with the central difference method. 

Comparison with Flight Test Data 

The longitudinal dynamic response is first 
calculated with the quasi-steady aerodynamic model, the 
pilot longitudinal control (S^ot) being the only input, in 
addition to the contributions from flight control systems. 

The net control deflections are plotted as 5^ 5„ and 5 r . The 
initial flight state variables are found to be the same as the 
initial flight data with a center of gravity location at 
x^O.45923 c . Figure 4 shows the simulation results in 
comparison with flight data. Although the longitudinal 
responses (a, 0 and q) follow the trend of flight data well, 
there are time periods in which the predicted a is much 
lower than the data. These are also the time periods with 
lower magnitude in pitch damping. The lateral-directional 
motions are not excited at all, although there are some 
initial inputs from the control system. The aforementioned 
longitudinal responses are repeated with the unsteady 
aerodynamic models as shown in Figure 5. Assuming that 
the control system in the simulation code models accurately 
the one in the flight vehicle, the discrepancy' in the 
responses must be caused by the static aerodynamic data 
which are also used in the unsteady models. The biggest 
differences between these two aerodynamic models occur in 
the lateral-directional responses. However, the excited 
lateral-directional motions with the unsteady models appear 
to be more unstable than what the flight data show, in 
particular in roll and sideslip 

Sideslipping, Roll 

To help explain the aerodynamic causes of some 
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loss-of-conlrol flight conditions, such as the USAir Flight 
427 accident in 1994°, the initial flight scenarios with a 
hard-over rudder malfunction will be postulated with the 
present F-16XL model. The present configuration is 
trimmed at M=0.1718 and 6000 ft. of altitude. With the 
longitudinal stick fixed at the trim position, the lateral- 
directional control systems are disconnected, but the 
longitudinal system is left on. In the initial 1 .5 seconds, a 
left bank is created by applying aileron and at the same 
time the rudder is set at a nose-left 20 deg.. After 1.5 
seconds, right aileron is applied to try to make the wing 
level. These are the scenarios indicated in Ref. 6. The 
results of simulation are presented in Figure 6. It is seen 
that with the quasi-steady aerodynamic model, a positive 
sideslip up to 12 deg. can be quickly generated while the 
angle of attack is increased to 36 deg. by the control system 
to counteract perhaps the increasing negative pitch attitude. 
Because of the positive sideslip, the resulting rolling 
moment overpowers the right aileron and the configuration 
keeps rolling to the left. Therefore, rapid sideslipping roll 
may generate loss-of-control flight conditions. 

If the unsteady aerodynamic models are used in 
the simulation, the situation gets worse, with the aircraft 
responses in a and P being divergent within five seconds. 
Typically, unsteady aerodynamics will make an unstable or 
marginally stable configuration with quasi-steady 
aerodynamic models more unstable. In the present case, the 
main reasons for the divergence are C ma being positive and 
C„p being negative as shown in Figure 7. Based on the 
present simulation, even maintaining a small rudder 
deflection will make the model highly unstable (not 
shown). Whether these highly unstable conditions with the 
unsteady aerodynamic models exist can only be verified 
with more appropriate test data. These needed data should 
be those obtained with coupled motions in pitch, roll and 
yaw. From a theoretical point of view, this is because the 
lateral-directional aerodynamics, such as roll and yaw 
damping, depend on wing boundaiy layer characteristics 
and the associated wake, which in turn strongly depend on 
the angle of attack and its time rate of change. Another 
implication is that these appropriate unsteady aerodynamic 
data may be needed in control system design 

Concluding Remarks 

Longitudinal and lateral-directional unsteady 
aerodynamic models obtained in separate wind tunnel 
testing were incorporated into an F-16XL simulation code 
which was based on a conventional quasi-steads 


aerodynamic model The resulting code was employed to 
investigate the differences in flight characteristics with 
these two types of aerodynamic models. For the unsteady 
model, the reduced frequency needed in the model was 
calculated with a least-square method to fit a harmonic 
variation to prior time history of state variables over about 
0.7 second. Methods to extract dynamic derivatives while 
performing time integration of dynamic equations were 
developed based on the concept of central differences and 
small-amplitude harmonic oscillation. These two methods 
showed similar results which were significantly different 
from the original quasi-steady data. Flight data in a 
longitudinal maneuver were compared with the simulation 
results. The main differences in dynamic response based on 
these two types of aerodynamic models were the lateral- 
directional responses induced by the longitu dinal motion. 
For the quasi-steady model, the lateral-directional motion 
was not excited, as compared with the flight data. On the 
other hand, the unsteady model produced a more unstable 
motion than what the flight data indicated. In a simulated 
sideslipping roll, again the configuration was shown to be 
more directionally unstable with an unsteady model. To 
improve the unsteady aerodynamic modeling, test data with 
coupled motions in pitch, roll and yaw would be needed. 
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Figure 1 Comparison between 
predicted results and assumed 
motion: 

a=35+20*cos(p/2*t) 
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Figure 2 Static aerodynamics 
and calculated dynamic derivatives 
based on samll amplitude harmonic 
oscillation(amplitude: a=5 deg) 
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Figure 3 Simulation results with unsteady aerodynamic model and quasi-steady 
aerodynamic model(pilot elevator input only). xcg=0. 45923 c 
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Figure 4 Comparison between flight test data and simulation data with 
quasi-steady aerodynamics at M=0.4243, H=14931.03ft and xcg=0 45923c 

10 

American Institute of Aeronautics and Astronautics 



§e.pilot( de 9) C Lq §e( de 9) 



time(sec) 



time(sec) 

10 
5 
0 

.5 t~i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 

0 5 10 15 20 25 30 35 40 45 

time(sec) 



time(sec) 




time(sec) 



time(sec) 



time(sec) 


Figure 4 Concluded 


11 


American Institute of Aeronautics and Astronautics 



fldeg) r(deg/s) p(deg/s) a(deg) 



time(sec) 



time(sec) 




time(sec) 



time(sec) 



time(sec) 


time(sec) 




time(sec) time(sec) 

Figure 5 Comparison between flight test data and simulation data with unsteady 
aerodynamic effects at M=0.4243, H= 1 4931 .03ft and xcg=0. 45923c 
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Figure 5 Concluded 
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Figure 6 Comparison between the simulation results with and without unsteady 
aerodynamic effects at M=0.1718, H=6000ft and xcg at 0.45923c 
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Figure 7 Stability derivatives of the unsteady configuration in flight 
conditions of Figure 6 


15 


American Institute of Aeronautics and Astronautics 



